#!/usr/bin/perl
# 2012/11/16 To use MACS, convert format to BED format using only potential self ligation product

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 4 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i  [map file] -o [output bed file]\n";
}

#---------------------------------------
# map format
#---------------------------------------
# id chr1 loc1 direction1 quality1 hind1 hind_loc1 chr2 loc2 direction2 quality2 hind2 hind_loc2 type(unique or repeat) number(alignable locations) locations

my %opt;
getopts("i:o:", \%opt);
my $FILE_map = $opt{i};	# include repeat, but no PCR bias data. Before removing self ligation products
my $FILE_out = $opt{o};

#---------------------------------------
# collect data
#---------------------------------------
my $SUM = 0;
my $NUM = 0;
my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";
my $fh_out = IO::File->new($FILE_out, 'w') or die "cannot write $FILE_out: $!";
while($_ = $fh_map->getline()){
	s/\r?\n//;
	my ($id, $chr1, $loc1, $direction1, $mapQ1, $hinID1, $hinLoc1, $chr2, $loc2, $direction2, $mapQ2, $hinID2, $hinLoc2, $type, $num, $locations) = split /\t/;


	# R or U
	my ($typeL, $typeR) = split /:/, $type;
	if($typeL ne 'U' and $typeL ne 'R'){
		next;
	}
	if($typeR ne 'U' and $typeR ne 'R'){
		next;
	}


	#---------------------------------------
	# check self ligation pattern
	#---------------------------------------
	my @lefts;
	push @lefts, [$chr1, $loc1,$direction1];
	my @rights;
	push @rights, [$chr2, $loc2, $direction2];

	my ($reL, $reR) = split /:/, $locations;
	if($typeL eq 'R'){
		if($reL eq 'NA'){
			next;
		}
		my @data = split /\;/, $reL;
		foreach my $k(@data){
			my ($chr, $L, $CIGR, $miss) = split /,/, $k;
			my $place = substr($L, 1);
			my $di = substr($L, 0, 1);
			push @lefts, [$chr, $place, $di];
		}
	}
	if($typeR eq 'R'){
		if($reR eq 'NA'){
			next;
		}
		my @data = split /\;/, $reR;
		foreach my $k(@data){
			my ($chr, $L, $CIGR, $miss) = split /,/, $k;
			my $place = substr($L, 1);
			my $di = substr($L, 0, 1);
			push @rights, [$chr, $place, $di];
		}
	}


	my $total_self_candidate = 0;
	foreach my $le(@lefts){
		my ($chr1_c, $loc1_c, $direction1_c) = @$le;
		foreach my $ri(@rights){
			my ($chr2_c, $loc2_c, $direction2_c) = @$ri;
			if(&checkSelfLigation($chr1_c, $chr2_c, $loc1_c, $loc2_c, $direction1_c, $direction2_c) == 1){
				$total_self_candidate++;
				$chr1 = $chr1_c;
				$chr2 = $chr2_c;
				$loc1 = $loc1_c;
				$loc2 = $loc2_c;
				$direction1 = $direction1_c;
				$direction2 = $direction2_c;
			}
		}
	}
	if($total_self_candidate != 1){
		next;
	}

	### please be careful about position
	my ($start, $end);
	if($direction1 eq "+"){
		$start = $loc2;
		$end = $loc1;
	}else{
		$start = $loc1;
		$end = $loc2;
	}

	$SUM += $end - $start;
	$NUM++;

	$fh_out->print("$chr1\t$start\t$end\t$id\t1\t$direction1\n");

}
$fh_map->close();
$fh_out->close();

my $average_fragment_size = $SUM / $NUM;
printf "Average fragment length: %.3f\n", $average_fragment_size;
printf "Shift size should be %d\n\n", int($average_fragment_size / 2);


sub checkSelfLigation{
	my ($chr1, $chr2, $loc1, $loc2, $direction1, $direction2) = @_;
	if($chr1 ne $chr2){
		return 0;
	}
	if(($direction1 eq '+') and ($direction2 eq '+') and ($loc2 < $loc1) and ($loc1 - $loc2 < 3000)){
		return 1;
	}
	if(($direction1 eq '-') and ($direction2 eq '-') and ($loc1 < $loc2) and ($loc2 - $loc1 < 3000)){
		return 1;
	}
	return 0
}


